 function [GG,impact,C,eu,SDX,ZZ,initss,ssvec,flag,ssnames,stanames,shonames]=modchimodDtFacpInfExpGapPistar(parest,options,op,calval,posest,poscal)
% =========================================================================
%
% Last Modified July 25 2011. Compatibility with code ion 
%
% function [GG,C,impact,eu,SDX,ZZ,NY,NX,NETA,gev,ncof]=...
% modchimodDtFacpXSpreadNews(parest,calval,posest,poscal)
%
%% *I. Description for CHIMOD DT FACP X SPREAD NEWS4* 
% 
% Model code for Chicago Fed model (CHIMOD) with 
% Double Trend (DT) Factor structure for prices (FACP) 
% Separates the ARMA(1,1) wage markup/labor disutility shock into a 
% white noise wage markup and an AR(1) labor disutility shock (X) 
% with a spread observable (SPREAD) and allowing for policy news shocks
% upto 4 periods ahead (NEWS4). 
% 
%% *II. Input* 
%
% There are either 1 or 4 inputs 
% *CASE 1* If there is 1 input, then PAREST should be the whole parameter
% vector. 
% 
% *CASE 2* If there are 4 inputs, the parameter vector PARAM is created
% from the estimated and calibrated values. 
% PAREST are the estimated values in rows POSEST and CALVAL are the
% calibrated values in position POSCAL. 
%
%% *III. Output*  
%
% The model to be solved is written, inside this code, as $$ A_{0}Y_{t+1} +
% A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  where the A matrices are
% filled with the linearized equation and entries depend on the parameter
% vector PARAM. 
%
% The solution of the model is a state space system characterized by 
%
% *Transition Equation* $$ S_{t} = GG S_{t-1} + RR \eta_{t} $$ and variance covariance matrix, 
% where $$ S_{t} $$ is the vector of all
% variables (endogenous and exogenous) of the model of dimension [NY 1],
% and $$ \eta $$  is the vector of exogenous innovations of dimension [NX 1] with 
% variance covariance matrix $$ V(\eta) = SDX'SDX $$. 
% Matrix dimensions: 
% GG [NY NY], IMPACT [NY NX], SDX [NX NX].  
%
% *Measurement Equation* $$ Y_{t} = ZZ S_{t} + ZZ C  $$ , 
% where $$ Y_{t} $$ is the data vector of dimension [NZ 1], $$ ZZ $$ is a selection matrix indicating 
% which states are observed, and $$ CC $$ is a vector of constants of
% dimension [NY 1]. 
% 
% EU is a [2 1] vector equal to [1;1] if there is a unique stable RE
% solution 
% 
% NETA GEV and COF are obsolete outputs that are still called in some codes
% and will be deleted in the new streamlined version of the estimation
% routine. 
%
%% *IV Notes on Multiple Prices* 
% 
% This version accomodates as observables an Investment Deflator and an additional NSERP prices. 
% Of these NSERP prices the first one should be the GDP deflator and the
% remaining NSERP-1 series multiple observables for the Consumption
% Deflator (e.g. chain weighted model consistent, Core PCE, Core CPI). 
% 
%% *V. Notes on Parametrization* 
% 1. Taylor rule written in terms of annualized output and inflation
% growth, with an inflation drift. 
% 2. Wage markup and labor disutility shock split. 
% 3. The composite steady state (SS) growth rate of the trend and the (SS)
% growth rate of ISTS progress are controlled by separate parameterss
% *gamma100* and *gammamiu100* rather than two parameters connected with
% ALPHA. This facilitates putting priors. 
% 4. AR coefficients of idiosyncratic errors are given by 2*PAR-1, where
% PAR is the paramater plugged into solution. Hence while PAR should be in
% the [0,1] interval, the AR coefficients are allowed to go negative. 
% 
%% *IV.a Model and data means* 
% 
% The model based means in the observation equation are given by the
% parameters *gammamiu100* and *pss100* which correspond to the mean growth
% rate of the relative price of investment to consumption and the mean
% growth rate of the consumption deflator. In addition, we allow the means
% of the observed price series to differ from these model based means
% according to extra parameters, intended to capture the discrepancy
% between model and data. 
% 
% Hence the model based mean of the Investment deflator is obtained as
% *pss100-gammamiu100-CMIU*.
% Similarly, the constants for the Consumption prices  is given by
% *pss100+CP* where CP differs by price. For simplicity this is also the
% constant for the GDP deflator, rather than aggregating that for
% investment and consumption; historically for instance Core PCE and the
% GDP deflator have grown at very similar rates. 
% 
%% *V. Related models* 
% 
% *MOD CHIMOD DT FAC P3 X SPREAD NEWS4 _ OBS:* is used to declare
% additional variables not used in estimation but annalized (e.g. output
% gap, real rate). 
% 
% *MOD CHIMOD DT FAC P3 SPREAD NEWS4:* very similar code but without split
% of wage markup and labor distulity shocks 
% 
%% *Version*
% Created by Alejandro Justiniano 12/13/2010  
% 
% Last Modified 12/13/2010. Run code ZB_CHECKWAGESPC to ensure alternative
% form of the wage Phillips curve matches model without split of wage markup and labor
% disutility shock. Also IRF of I to neutral shock is reasonable. 
% =========================================================================

if nargout < 11 
    flag_repstanames=0; 
else 
    flag_repstanames=1; 
end 
if nargout < 10 
    flag_repssnames=0; 
else
    flag_repssnames=1; 
end

flag.ssok=1;
flag.solver=1; 

initss=[]; 
ssnames=[]; 
ssvec  =[]; 


% if nargin~=1 && nargin~=4 
%     error('Need either 1 or 4 inputs') 
% end 

% =========================================================================
%% *PART 1. Declare All Variables in the Model Solution $$ S_{t} $$* 

%__________________________________________________________________________
% 1.a Endogenous Variables 
y=1;            % output
k=2;            % capital
L=3;            % hours
Rk=4;           % return on capital
w=5;            % real wages
x=6;            % marginal utility of labor
p=7;            % inflation
s=8;            % marginal cost
lambda=9;       % multiplier
c=10;           % consumption
R=11;           % interst rate
u=12;           % capital utilization
phi=13;         % multiplier
i=14;           % investment
kbar=15;        % "gross" capital
%__________________________________________________________________________
% 1.b Exogenous driving processes
mp=16;
z =17;          % Productivity shock
g =18;          % Public spending shock
miu=19;         % Relative price (non stationary), this is Growth Rate
lambdap=20;     % Markup shock AR(1) 
psi=21;         % Preference shock AR(1) 
b=22;           % Discount factor shock 
upsilon=23;     % MEI, stationary, investment shock
pitarg=24;      % Inflation Drift
wmarkwn=25;     % White Noise Wage Markup 
%__________________________________________________________________________
% 1.c Observable growth rates (nominal) 
ynomobs=26; 
cnomobs=27;
inomobs=28; 
wnomobs=29; 
%__________________________________________________________________________
% 1.c.a Auxiliary Variables or Observable Prices
yrobs   =30; 
gdp     =31; 
dpinv   =32; 
gdpdef  =33; 
%__________________________________________________________________________
% 1.d Lags
yrobs_1=34;   % Real Observable Output (adjusted for utilization) lag1      
yrobs_2=35;   % Real Observable Output (adjusted for utilization) lag2
yobsan =36;   % Annualize Observable Output 
p_1    =37; 
p_2    =38; 
dpma   =39;   % Four quarter inflation  
dagtrend=40;  % First difference of the aggregate trend 
%__________________________________________________________________________
% 1.e 1 through 4 step ahead expectations of FFR 
startstar=40;
% 1.f Policy signals
stsig=startstar; 
sigR_0=stsig+1; 
sigR_1=stsig+2; 
sigR_2=stsig+3; 
sigR_3=stsig+4; 
sigR_4=stsig+5; 
% 1.g Inflation Expectations 
% posExpInf=stsig+5+1; 
% %% Expectations of inflation 1 through 40 quarters ahead
% pe_vec=(posExpInf):(posExpInf+39); 
% expectedPi40=posExpInf+40; 
NY=stsig+5;

% =========================================================================

% =========================================================================
%% *PART 2. Declare Innovations to the exogenous variables $$ \eta_{t} $$* 

nsignals = 4;    % Number of policy signals 
NX=10+ nsignals; 
%__________________________________________________________________________
% 2.a Structural Disturbances 
Rs=1;           % MP
zs=2;           % Technology 
gs=3;           % Public spending
mius=4;         % Investment     
lambdaps=5;     % Mark-up
psis=6;         % Leisure preference
bs=7;           % Discount factor
upsilons=8;     % investment shock
shpitarg=9;     % Inflation Drify
wmarkwnsh=10; 
%__________________________________________________________________________
% 2.b Policy Signals 
Rs_lag1=11;    % 1 period ahead
Rs_lag2=12;    % 2 periods ahead
Rs_lag3=13;    % 3 periods ahead
Rs_lag4=14;    % 4 periods ahead
% =========================================================================

% =========================================================================
%% *PART 3. Declare Parameters* 
ncof=        69;    % Number of coefficients 
numpar=      ncof;  % Equal to Parameters, distinction irrelevant here 
if nargin < 4 || isempty( calval ) 
    param=parest  ; 
else
    param=zeros(numpar,1); 
    param(posest)=parest; 
    param(poscal)=calval; 
end; 
%__________________________________________________________________________
% 3.1 Parameters for all but multiple prices, GDP and I deflators 

alpha=param(1);         % share of capital in the prod. function
delta=param(2);         % capital depreciation rate
iotap=param(3);         % price indexation (=0 is static indexation, =1 is dynamic)
iotaw=param(4);         % wages indexation
gammastar100=param(5);  % steady state growth rate of technology
gammamiu100=param(6);   % steady state growth rate of capital embodied technology
h=param(7);             % habit formation parameter
lambdapss=param(8);     % steady state of mark-up shock (pins down steady state of wages)
lambdaw=param(9);       % wage mark-up
Lss=param(10);          % steady state for leisure (the ss for leisure is pinned down by psiss. convenient to parameterize the model in terms of Lss)
pss100=param(11);       % steady state quarterly inflation (multiplied by 100)
Fbeta=param(12);        % weird parameter of SW that is a function of beta and the steady state quarterly real rate of interest (ensures that beta is less than 1)
gss=param(13);          % steady state of the public spending shock
niu=param(14);          % curvature of the utility function for leisure
xip=param(15);          % price stickiness
xiw=param(16);          % wage stickiness
chi=param(17);          % elasticity of the capital utilization cost function
S=param(18);            % investment adjustment cost
fp=param(19);           % reaction to annualized inflation
fy=param(20);           % reaction to annualized output growth 
rhoR=param(21);         % policy smoothing  
% AR coefficients of exogenous prorcesses 
rhoz=param(22); 
rhog=param(23); 
rhomiu=param(24); 
rholambdap=param(25); 
rhopsi=param(26); 
rhob=param(27); 
rhopitarg=param(28); 
rhomp=param(29); 
rhoupsilon=param(30); 
% Std of structural shocks 
sdstruct=param(31:40);
kappa       = -param(41);   % Scaling for spredd  
sdspreadme  =  param(42);   % Measurement error spread
% Std Policy Signals
stdsignals  =param(43:46); 
if length(stdsignals)~=nsignals 
    error('Wrong STD for signals') 
end 
%% Parameters for inflation expectations 
posCPI=param(47); 
expInfLoading=param(48); 
expInfCons   =param(49); 
end_nonp=49; 
% _________________________________________________________________________
% 3.2  Begin Parameters for Multiple Prices 
%__________________________________________________________________________
% 3.2.1  How to aggregate gdp deflator
flag_wadd1  =param(end_nonp+1); 
nserp       =param(end_nonp+2); 
% __________________________________________________________________
% 3.2.1 Declare positions for 
% Rows for 
% a) loadings 
% b) ar ID errors 
% c) constants 
% d) volatilities 
%
% *Note:* There are NSERP+1 loadings. First 2 are for the GDP Deflator.
% 1st of W(1)*Price consumption.2nd on W(2)*Prince Investment. 
% Remaining are for other prices on the Consumption Deflators 
row_load=(end_nonp+3):(end_nonp+nserp+3); 
row_ar  =row_load(end)*ones(1,nserp)+(1:nserp); 
row_std =row_ar(end)*ones(1,nserp)+(1:nserp);
row_c   =row_std(end)*ones(1,nserp)+(1:nserp); 

%__________________________________________________________________________
% 3.2.3 vector of factor loadings 
bet_gdpdef  =param(row_load(1:2)); 
bet_c       =param(row_load(3:end)); 

bet_CPI     =bet_c(posCPI-1); 
if abs(expInfLoading-999) < 0.01; 
    expInfLoading=bet_CPI;
end



%__________________________________________________________________________
% 3.2.4 vector of idiosyncratic errors 
rhoid=2*param(row_ar)-1; 
if any(rhoid < -0.995)~=0 || any(rhoid > 0.995)~= 0 
    disp('Explosive AR ID roots'); 
    GG=[];C=[];impact=[];eu=[0;0];SDX=[];ZZ=[];NY=[];NX=[];NETA=[];gev=[];ncof=[];
    return
end

%__________________________________________________________________________
% 3.2.5 vector of constants 
cprices =param(row_c); 
if all( cprices ~= 0 ) 
    error('At least one of the constants must be set equal to 0')
end
%constantCPI=cprices(posCPI); 

%__________________________________________________________________________
% 3.2.6 vector of idiosyncratic errors
sdid             =param(row_std); 

%__________________________________________________________________________
% 3.2.7 constants for I deflator 
cons_idef        =param(row_c(end)+1:end);
if length(cons_idef)> 1 
    error('Wrong Dimension of Investment Deflator') 
end 

%==========================================================================
%% *PART 4. Compute Steady State (SS)* 
 
gammastar=gammastar100/100;
gammamiu=gammamiu100/100;
beta=100/(Fbeta+100);
rss=exp(gammastar)/beta-1;
rss100=rss*100;
gss=1/(1-gss);
Rkss=exp(gammastar+gammamiu)/beta-1+delta;
sss=1/(1+lambdapss);
wss=(sss*((1-alpha)^(1-alpha))/((alpha^(-alpha))*Rkss^alpha))^(1/(1-alpha));
kLss=(wss/Rkss)*alpha/(1-alpha);
FLss=(kLss^alpha-Rkss*kLss-wss);
yLss=kLss^alpha-FLss;

% Scale of the economy is fixed at 1, i.e. independent of LSS which now
% only enters observation equation for hours 
Lsscale=1;
kss=kLss*Lsscale;
iss=(1-(1-delta)*exp(-gammastar-gammamiu))*kss*exp(gammastar+gammamiu);
F=FLss*Lsscale;
yss=yLss*Lsscale;
css=yss/gss-iss;

ssvec  =[css/yss;iss/yss;kss/yss;rss100;100*Rkss]; 
if  flag_repssnames==1;
    ssnames={'c/y'  ;'i/y'  ;'k/y';  'rss100';'Rkss100'};
end
    
%%==========================================================================

%==========================================================================
%% *PART 5. Declare System Matrices* 
% Model to be solved is written as 
% $$ A_{0}Y_{t+1} + A_{1}Y_{t}=A_{2}Y_{t-1} + PSI \eta_{t} $$  
A0  = zeros(NY,NY) ;
A1  = zeros(NY,NY) ;
A2  = zeros(NY,NY) ;
PSI = zeros(NY,NX) ;
%==========================================================================
expg=exp(gammastar);

%% *PART 6. Begin Filling Equation Rows* 

%% Eq 1, Output (y) from production function 
A1(y,y)=1;
A1(y,k)=-((yss+F)/yss)*alpha;
A1(y,L)=-((yss+F)/yss)*(1-alpha);

%% Eq. 2, Effective Capital (K) 
A1(k,k)=1;
A1(k,u)=-1;
A1(k,z)=1;
A2(k,kbar)=1;
A1(k,miu)=alpha/(1-alpha)+1;

%% Eq. 3 Hours (L), from cost minimization 
A1(L,Rk)=1;
A1(L,w)=-1;
A1(L,k)=1;
A1(L,L)=-1;

%% Eq. 4 Return to Capital (Rk) from utilization FOC
A1(Rk,Rk)=1;
A1(Rk,u)=-chi;

%% Eq. 5, Real Wages (w) PC (with split of shocks)
slopewages=(1-beta*xiw)*(1-xiw)/((1+beta)*xiw*(1+niu*(1+1/lambdaw)));
A1(w,w)=1;
A0(w,w)= -beta/(1+beta);
A0(w,p)= -beta/(1+beta);
A1(w,x)=slopewages;
A1(w,z)=(1+iotaw*beta-rhoz*beta)/(1+beta);
A1(w,miu) =(alpha/(1-alpha))*((1+iotaw*beta-rhomiu*beta)) *(1/(1+beta)); 
A1(w,p)=(1+iotaw*beta)/(1+beta);
A1(w,wmarkwn)= -1; 

A2(w,p)=iotaw/(1+beta);
A2(w,w)=1/(1+beta);
A2(w,z)=iotaw/(1+beta);
A2(w,miu)=(alpha*iotaw)/( (1+beta)*(1-alpha) ); 

%% Eq. 6, Wage Gap (x), (MRS-Real Wage)
A1(x,x)=1;
A1(x,w)=-1;
A1(x,lambda)=-1;
A1(x,b)=1/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2));
A1(x,L)=niu;
A1(x,psi)=1; 

%% Eq 7, Prices (p) PC 
A1(p,p)=1;
A0(p,p)=-beta/(1+iotap*beta);
A2(p,p)=iotap/(1+iotap*beta);
A1(p,s)=-((1-beta*xip)*(1-xip)/((1+iotap*beta)*xip));
A1(p,lambdap)=-1;       

%% Eq 8, Marginal Cost (s)
A1(s,s)=1;
A1(s,Rk)=-alpha;
A1(s,w)=-(1-alpha);

%% Eq 9, Lagrange Multiplier IBC (Lambda) 
A1(lambda,lambda)=1;
A1(lambda,R)=-1;
A0(lambda,lambda)=-1;
A0(lambda,p)=1;
A1(lambda,z)=rhoz;
A1(lambda,miu)=(rhomiu)*alpha/(1-alpha);

%% Eq 10, Consumption (c)
A1(c,lambda)=(expg-h*beta)*(expg-h);
A1(c,b)=-(expg-h*beta*rhob)*(expg-h)/((1-rhob)*(expg-h*beta*rhob)*(expg-h)/(expg*h+expg^2+beta*h^2)); 
A1(c,z)=-(beta*h*expg*rhoz-h*expg);
A1(c,miu)=(-(beta*h*expg*rhomiu-h*expg))*alpha/(1-alpha);
A1(c,c)=(expg^2+beta*h^2);
A0(c,c)=-beta*h*expg;
A2(c,c)=expg*h;

%% Eq 11, Taylor Rule for Nominal Interest Rate (R)
A1(R,R)=1;
% 10.1 Systematic Component 
A2(R,R)=rhoR;
A1(R,dpma)  = -(1-rhoR)*fp; % Annualized Inflation 
A1(R,gdp )  = -(1-rhoR)*fy; % Current Gap
% 10.2 Non-Systematic Component 
A1(R,mp)      = -1;         % Contemporaneous Innovation 
A1(R,sigR_0)  = -1;           % Signals 
A1(R,pitarg)  = expInfLoading; % Inflation Drift

%% Eq. 12, Utilization (U) from market clearing condition 
A1(u,c)=css/yss;
A1(u,i)=iss/yss;
A1(u,y)=-1/gss;
A1(u,g)=1/gss;
A1(u,u)=kss*Rkss/yss;

%% Eq. 13, Shadow Value of Capital (phi)
A1(phi,phi)=1;
A0(phi,phi)=-beta*exp(-gammastar-gammamiu)*(1-delta);
A1(phi,z)=rhoz;
A1(phi,miu)=(rhomiu)*(alpha/(1-alpha)+1);
A0(phi,lambda)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));
A0(phi,Rk)=-(1-beta*exp(-gammastar-gammamiu)*(1-delta));

%% Eq. 14, Investment (i, in consumption units) 
expgmiu=exp(gammastar+gammamiu);
A1(i,lambda)=1/(S*expgmiu^2);
A1(i,phi)=-1/(S*expgmiu^2);
A1(i,upsilon)=-1/(S*expgmiu^2);
A1(i,miu)=((1-beta*rhomiu))*(alpha/(1-alpha)+1); %no normalization
A1(i,i)=(1+beta);
A1(i,z)=(1-beta*rhoz);
A0(i,i)=-beta;
A2(i,i)=1;

%% Eq 15, Capital accumulation (kbar)
% Note: Recall, effective capital above is KBAR*UTILIZATION 
A1(kbar,kbar)=1;
A1(kbar,i)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A1(kbar,upsilon)=-(1-(1-delta)*exp(-gammastar-gammamiu));
A2(kbar,kbar)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,z)=(1-delta)*exp(-gammastar-gammamiu);
A1(kbar,miu)=((1-delta)*exp(-gammastar-gammamiu))*(alpha/(1-alpha)+1);

%% Eq 16, MP Exogenous Policy Innovation (unexpected) 
A1(mp,mp)=1; A2(mp,mp)=rhomp; PSI(mp,Rs)=1; 

%% Eq 17, Exogenous Growth Rate in Neutral Technoloy 
A1(z,z)=1; A2(z,z)=rhoz; PSI(z,zs)=1;

%% Eq 18, Exogenous Process in G Spending 
A1(g,g)=1; A2(g,g)=rhog; PSI(g,gs)=1;

%% Eq 19, Exogenous Process in Investment Specific Tech or Relative Price C
% to I 
A1(miu,miu)=1; A2(miu,miu)=rhomiu; PSI(miu,mius)=1;

%% Eq 20, Exogenous Process in Price Markup 
A1(lambdap,lambdap)=1; A2(lambdap,lambdap)=rholambdap;PSI(lambdap,lambdaps)=1;

%% Eq 21, Exogenous Process in Labor Disutility (AR1) 
A1(psi,psi)=1; A2(psi,psi)=rhopsi; PSI(psi,psis)=1;

%% Eq 22, Exogenous Process in Intertemporal Preference
A1(b,b)=1; A2(b,b)=rhob; PSI(b,bs)=1;

%% Eq 23, Exogenous Process in MEI 
A1(upsilon,upsilon)=1; A2(upsilon,upsilon)=rhoupsilon; PSI(upsilon,upsilons)=1;

%% Eq 24, Exogenous Process in Inflation Drift 
A1(pitarg,pitarg)=1; A2(pitarg,pitarg)=rhopitarg; PSI(pitarg,shpitarg)=1;

%% Eq 25, Exogenous Process in Wage Markup (white noise) 
A1(wmarkwn,wmarkwn)=1; PSI(wmarkwn,wmarkwnsh)=1; 

%% Eq. 26, Observable GDP Nominal (add P back)
A1(ynomobs,ynomobs) =  1;
A1(ynomobs,gdp)     = -1; 
A1(ynomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(ynomobs,p)       = -1; 

A2(ynomobs,gdp)     = -1; 

%% Eq. 27, Observable Consumption Nominal (add P back)
A1(cnomobs,cnomobs) =1;
A1(cnomobs,c)       = -1; 
A1(cnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(cnomobs,p)= -1; 

A2(cnomobs,c)= -1; 

%% Eq. 28 Observable Investment Nominal (add P back)
% In consumption units
A1(inomobs,inomobs)   =1;
A1(inomobs,i)         = -1; 
A1(inomobs,dagtrend)  = -1; 
% Add inflation, since nominal 
A1(inomobs,p)= -1; 

A2(inomobs,i)= -1; 

%% Eq. 29 Observable Hourly Wages  Nominal (add P back)
A1(wnomobs,wnomobs) =1;
A1(wnomobs,w)       = -1; 
A1(wnomobs,dagtrend)= -1; 
% Add inflation, since nominal 
A1(wnomobs,p)= -1; 

A2(wnomobs,w)= -1; 

%% Eq. 30 Observable Real GDP Growth (in C Units)
% Corrected observation equation for utilization 
% $$ yrobs_{t} = gdp_{t}-gdp_{t-1} + agtrend_{t} $$ 
% where AGTREND is the aggregate trend 
A1(yrobs,yrobs)=1; 
A1(yrobs,gdp)= -1;
A1(yrobs,dagtrend)= -1; 
A2(yrobs,gdp)= -1; 

%% Eq. 31 GDP (detrended), Real Detrended Output corrected for Utilization 
A1(gdp,gdp)= -1; 
A1(gdp,y)=1; 
A1(gdp,u)= -kss*Rkss/yss; 

%% Eq. 32 Growth Rate in Investment Deflator 
% $$ \Delta P^{I}_{t} = \Delta P^{C}_{t} - \mu_{t} $$
A1(dpinv,dpinv)= -1; 
A1(dpinv,p)    =  1; 
A1(dpinv,miu)  = -1; 

%% Eq. 33 GDP Deflator 
% $$ P^{GDP}_{t} = \beta_{1} P^{C}_{t} + \beta_{2} P^{I}_{t} $$ 
% 
% Recall there are 2 ways to define the weights 
%
% *First:* $$ {\frac{C}{C+I} , \frac{I}{C+I} } $$ 
% 
% *Second:* $$ {\frac{C}{Y}  , \frac{I}{Y}   } $$
%
% Note that in second case weights will not add up to one, which is not
% desirable, but opens up the door for adding a 3 deflator, i.e. price of G
% and NX. 
bet_gdpdef=bet_gdpdef(:); 
weights=zeros(1,2); 
if flag_wadd1==1
    weights(1)=css/(css+iss);
    weights(2)=iss/(css+iss);
else
    weights(1)=css/yss;
    weights(2)=iss/yss;
end
A1(gdpdef,gdpdef) = -1; 
A1(gdpdef,p)      = weights(1)*bet_gdpdef(1); 
A1(gdpdef,dpinv)  = weights(2)*bet_gdpdef(2); 

%% Eq. 34 Lag 1 Observable Real GDP Growth (in C Units)
A1(yrobs_1,yrobs_1)=1; A2(yrobs_1,yrobs)=1;

%% Eq. 35 Lag 2 Observable Real GDP Growth (in C Units)
A1(yrobs_2,yrobs_2)=1; A2(yrobs_2,yrobs_1)=1;

%% Eq. 36 Annualized Observable Real GDP Growth (in C Units)
% $$ \Delta y^{A}_{t}= \sum_{j=0}^{3}yrobs_{t-j} $$
A1(yobsan,yobsan)= -1;

A1(yobsan,yrobs)  =1; 
A1(yobsan,yrobs_1)=1; 
A1(yobsan,yrobs_2)=1;
A2(yobsan,yrobs_2)= -1;

%% Eq. 37 Lag 1 Inflation (C)
A1(p_1,p_1)=1; A2(p_1,p)=1;

%% Eq. 38 Lag 2 Inflation (C)
A1(p_2,p_2)=1; A2(p_2,p_1)=1;

%% Eq. 39 Annualized Inflation (C) 
% $$ \Delta P^{A,C}_{t}= \frac{ \sum_{j=0}^{3} \Delta P^{C}_{t-j} }{4} $$
A1(dpma,dpma)= -1;
A1(dpma,p)   = 1/4;
A1(dpma,p_1) = 1/4;
A1(dpma,p_2) = 1/4;
A2(dpma,p_2) = -1/4;

%% Eq. 40, Growth Rate in Aggregate Trend 
% $$ \Gamma_{t} = Z_{t} + \frac{\alpha}{1-\alpha} \mu_{t} $$ 
A1(dagtrend,dagtrend)= -1; 
A1(dagtrend,z)=1; 
A1(dagtrend,miu)=alpha/(1-alpha); 

%% Eq. 45-49 Monetary Policy Signals

%% Eq. 45. $$ \zeta^{4}_{t}       = \epsilon^{4}_{t} $$ 
A1(sigR_4,sigR_4)=1;PSI(sigR_4,Rs_lag4)=1; 

%% Eq. 46.  $$ \zeta^{3}_{t} = \zeta^{4}_{t-1} + \epsilon^{3}_{t} $$ 
A1(sigR_3,sigR_3)=1; A2(sigR_3,sigR_4)=1; PSI(sigR_3,Rs_lag3)=1; 

%% Eq. 47.  $$ \zeta^{2}_{t} = \zeta^{3}_{t-1} + \epsilon^{2}_{t} $$ 
A1(sigR_2,sigR_2)=1; A2(sigR_2,sigR_3)=1; PSI(sigR_2,Rs_lag2)=1; 

%% Eq. 48.  $$ \zeta^{1}_{t} = \zeta^{2}_{t-1} + \epsilon^{1}_{t} $$ 
A1(sigR_1,sigR_1)=1; A2(sigR_1,sigR_2)=1; PSI(sigR_1,Rs_lag1)=1; 

%% Eq. 49. $$ \zeta^{0}_{t} = \zeta^{1}_{t-1} $$, which implies that the history of
% policy signals regarding the current quarter is given by 
%
% $$ \zeta^{0}_{t} = \epsilon^{1}_{t-1} +  \epsilon^{2}_{t-2} + \epsilon^{3}_{t-3} +
% \epsilon^{4}_{t-4}   $$ 
% *Note:* The contemporaneous policy innovation appears 
% directly in the definition of the Taylor rule, in case it is AR(1). An
% alternative would be to make the current composite signal +
% contemporaneous innovation persistent. 
A1(sigR_0,sigR_0)=1;A2(sigR_0,sigR_1)=1;  


%% 
% pe_vec=pe_vec(:); 
% A1(pe_vec,pe_vec)=eye(length(pe_vec)); 
% A0(pe_vec,[p;pe_vec(1:end-1)])= -eye(length(pe_vec)); 
% 
% A1(expectedPi40,expectedPi40)= -1; 
% A1(expectedPi40,pe_vec      )= (1/length(pe_vec))*ones(1,length(pe_vec));  
% 

%% *PART 7: Solve Model using AMASOLVE*
[GG,~,impact,~,~,~,~,eu]=amasolve(A0,A1,A2,zeros(NY,1),PSI) ;
% EU =[1;1] if unique and stable RE solution 
if ~isequal(eu,[1;1]); 
    flag.euok=0; 
    ZZ=[]; C=[]; SDX=[];
    return 
end 

flag.euok=1; 
%% *PART 8: Addition of Measurement Equations for Prices*

%% 8.1 Extend model solution by 2*NSERP rows 
% Number of rows to add to the solution 
nadd=nserp*2; 
% Part 1: Position of prices 
posp =(NY+1):(NY+nserp); 
% Part 2: Position of ID errors 
posid=(posp(end)+1):(NY+2*nserp); 
GG=blkdiag(GG,zeros(nadd)); 
% Matrix of impact coefficients 
RR=zeros(NY+nadd,NX+nserp); 
RR(1:NY,1:NX)=impact; 

%% 8.2 Idiosyncratic errors 
% Which evolve as $$ \nu^{i}_{t} = \rho  \nu^{i}_{t-1} + \varsigma^{i}_{t} $$
% as independent AR(1) processes 
GG(posid,posid   )=diag(rhoid); 
RR(posid,NX+1:end)=eye(nserp); 

%% 8.3 General Structure of Observable Price Equation 
% 
% For price $$ P^{i,C}_{t} = \beta^{i} P^{C}_{t} + \nu^{i}_{t} $$ written in terms
% of lagged variables as 
% $$ P^{i,C}_{t} = \beta^{i} G_{P^{C}} S_{t-1} + \beta^{i} R_{P^{C}}
% \eta_{t} + \rho \nu^{i}_{t-1} + \varsigma^{i}_{t} $$ ,
% where $$ \beta^{i} $$ is the series specific factor loading, and, $$ G_{P^{C}}  R_{P^{C}} $$ correspond to the rows of the solution
% matrix for the model consumption price 

%% 8.3.a First row is GDP 
GG(posp(1),1:NY)=GG(gdpdef,1:NY);
RR(posp(1),1:NX)=RR(gdpdef,1:NX); 

%% 8.3.b Remaining rows are multiple indicators for consumption prices 
% Recall Structural part of observation equation $$ \beta^{i} G_{P^{C}} S_{t-1} +
% \beta^{i} R_{P^{C}} \eta_{t}  $$ 
GG(posp(2:end),:)       =repmat( GG(p,:)    ,[nserp-1 1] ); 
RR(posp(2:end),1:NX)    =repmat( RR(p,1:NX) ,[nserp-1 1] );
% Multiply by Factor Loadings  
bet_c=bet_c(:); 
GG(posp(2:end),1:NY)=(repmat(bet_c,[1 NY])).*GG(posp(2:end),1:NY); 
RR(posp(2:end),1:NX)=(repmat(bet_c,[1 NX])).*RR(posp(2:end),1:NX);
%% 8.3.c Add idiosyncratic errors 
GG(posp,posid)   =diag(rhoid); 
RR(posp,NX+1:end)=eye(nserp); 

% Update model dimensions after adding prices 
impact=RR; 
NY=size(GG,1); 
NX=size(RR,2); 

%% *Part 9: Addition of Measurement Equations for Spread*
% 
%% 9.1 Extend model solution by 1 row 
GG    =blkdiag(GG,0);
impact=[impact zeros(NY,1);zeros(1,NX+1)];

%% 9.2 Fill in last row for the spread 
% 
% $$ sp_{t} = \kappa \upsilon_{t} + \varsigma^{sp}_{t} $$ 
posspread=NY+1;
GG(posspread,:)       =kappa*GG(upsilon,:); 
impact(posspread,:)   =kappa*impact(upsilon,:); 
impact(posspread,NX+1)=1; 

% Update model dimensions after adding spread
NY=NY+1;
NX=NX+1;

%% *PART 10: Declare Cholesky of Variance Covariance Matrix* 
% Recall $$ V = SDX'SDX $$ and the order of the shocks is 
% 1. Structural 2. Policy Signals 3. Idiosyncratic Prices 4. Measurement
% Error Spread 
SDX=diag([sdstruct(:);stdsignals(:);sdid(:);sdspreadme]);
if size(SDX,1)~=NX;error('Dimensions of errors and SDX do not match');end

%% *PART 11: Declare ZZ, Matrix of Pointers to Observables*
% Ordering 
% 1) Nominal GDP Growth. 
% 2) Nominal C   Growth.
% 3) Nominal I   Growth.
% 4) Hours. 
% 5) Nominal W   Growth. 
% 6) Feds Fund Rate. 
% 7) Investment Deflator. 
% 8) GDP Deflator. 
% 9 though NSERP-1) Indicators of C Prices. 
% Last is Spread
dimZ=7+nserp+1+1; 
ZZ=zeros(dimZ,NY); 
ZZ(1,ynomobs)=1; 
ZZ(2,cnomobs)=1; 
ZZ(3,inomobs)=1; 
ZZ(4,L)=1; 
ZZ(5,wnomobs)=1; 
ZZ(6,R)=1;
ZZ(7,dpinv)=1;
ZZ(8:dimZ-2,posp)=eye(nserp);
ZZ(dimZ-1,posspread)=1; 
ZZ(dimZ  ,pitarg)=1; 

%% *PART 12: Declare C, Matrix of Constants*
C=zeros(NY,1); 
C(ynomobs)=gammastar100+pss100;
C(cnomobs)=gammastar100+pss100;
C(inomobs)=gammastar100+pss100;
C(L)=Lss;
C(wnomobs)=gammastar100+pss100;
C(R)    =pss100+rss100;
C(dpinv)=(pss100-gammamiu100)+cons_idef; 
C(posp )= pss100*ones(nserp,1)+cprices(:); 
C(pitarg)=pss100+expInfCons;

if flag_repstanames==0; 
    return 
end 
%________________________________________________________
%% 16.1 STANAMES: Cell with Names of ALL states 
stanames1={'yhat','khat','hours','Return K','what','x','dp','Marginal Cost','lambdahat','chat',...
    'nomr','Utilization ','phihat','ihat','kbar','mp','z','g','miu','lambdap',...
    'psi','b','upsilon ','pitarg','wmarkwn','GDPm nominal','Cm nominal','Im nominal','Wm nominal',...  % Position 29
    'GDPm real','GDPm level','I deflator','GDPm deflator','GDPm real lag1','GDP real lag2','GDP Annualized','Inflation lag1','Inflation lag2','Inflation annualized','Aggregate trend',...
    'Signal0','Signal1','Signal2',...
    'Signal3','Signal4',...    
    'GDP deflator'}; 
if nserp==3 
    stanames2={'C deflator','PCE Core','GDP def Id','C def Id','PCE Core Id'}; 
elseif nserp==4 
    stanames2={'C deflator','PCE Core','CPI Core','GDP def Id','C def Id','PCE Core Id','CPI Core Id'}; 
else 
    disp('Nprices > 4 please add names to STANAMES') 
end 
%stanames3={'Spread','Spread Me','RealY','RealC','RealI','RealW','SumR'}; 
stanames=[stanames1(:);stanames2(:);'Spread']; 

%% 16.2 STSHOCKS: Cell with Names of Shocks in STANAMES 
% *Note:* Technically not true for the policy signals, for which looking
% only at the cummulative signals. Hence isolate policy innovations rather
% than policy states.
stshocks={'mp','z','g','miu','lambdap','psi','b','upsilon ','pitarg','wmarkwn',...
    'Signal1','Signal2','Signal3','Signal4','GDP def Id','C def Id','PCE Core Id'};
if nserp==4
    stshocks=[stshocks,'CPI Id'];
end
stshocks=[stshocks,'Spread Me'];
stshocks=stshocks(:);

%% 16.3 SHONAMES: Cell with Names of Shock Innovations 
shonames={'MP','Neutral','G+NX','ISTS','Price Markup','Labor Disutility',...
    'Discount','MEI','Inflation Drift','Wage Markup',...
    'MP Signal 1','MP Signal 2','MP Signal 3','MP Signal 4'}; 
shonames=[shonames(:);stshocks(15:end)]; 


%obsnames={}; for ii=1:size(ZZ,1); temp=stanames( find(ZZ(ii,:)~=0 ) ); obsnames=[obsnames;temp(:)]; end;

%% 16.4 Quick check dimensions 
if length(stanames)~=size(GG,1)
    error('Mistmatch STANAMES and NY');
end
if length(stshocks)~=size(SDX,1)
    error('Mistmatch STSHOCKS and NX');
end
if length(shonames)~=size(SDX,1)
    error('Mistmatch SHONAMES and NX');
end
